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We extend the high-temperature series of the free energy for the XY model in two 
dimensions to order /3 48 from the previous order of /3 22 by applying an improved algorithm 
of the finite lattice method. The long series obtained allows us to conclude that the behavior 
of the free energy is consistent to high accuracy with what is expected when the phase 
transition of the model is of the Kosterlitz-Thouless type. 

§1. Introduction 

We investigate the XY model in two dimensions. It is believed to exhibit a 
phase transition of the Kosterlitz-Thouless (K-T) typep^ driven by the condensation 
of vortices. Through use of renormalization group arguments, it was predicted that 
the correlation length possesses an essential singularity at the transition temperature 
T c of the form 

' b 1 

where t = T/T c — 1 is the reduced temperature and b is a non-universal constant. 
The free energy and its temperature derivatives (i.e., the internal energy and the 
specific heat) have also been predicted to behave as 



/ ~ exp 



2b 

7T/2 



+ [ a regular term]. (T2) 



The first term on the right-hand side of Eq. (jl-2p has an essential singularity at t = 0. 
This term itself and its derivatives are zero at the critical temperature. The second 
term is regular at the critical temperature. 

The behavior expressed in (|l-ip for the correlation length has been well es- 
tablished both by numerical simulations and the high-temperature expansion. The 
standard Monte Carlo simulation^ gives b = 2.15(10) with (3 C = 1.130(15) and 
b = 1.70(20) with (3 C = 1.118(5) for the square lattice. (Here [5 C is the inverse critical 
temperature, which is defined below.) The more precise values b = 1.800(2) with 
P c = 1.1208(2) and b = 1.776(4) with (3 C = 1.1199(1) were obtained using the finite- 
size scaling technique® and the renormalization group finite-size scaling method,^ 
respectively. In the latter approach, the renormalization group flow of the observable 
was matched with that of the exactly solvable BCSOS model. The latter value of 
(3 C was recently confirmed by large scale Monte Carlo simulations on a 2048 x 2048 
lattice using the finite-size scaling methodP 1 A high-temperature expansion for the 
correlation lengtrP® gives the slightly smaller value b = 1.67(4) with f3 c = 1.118(3). 
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In contrast to the situation for (jl-lh . the behavior expressed by (11-21) for the 
free energy or its temperature derivatives has not yet been confirmed in numerical 
simulations or in high-temperature expansions. The reason for this is the following. 
The free energy and its derivatives are dominated by the regular term near the 
critical temperature, and the signal of the singularity in the first term is relatively 
weak. For this reason, a much more precise simulation is needed to confirm the 
behavior given in (jl-2f) for these quantities than in the case of the correlation length. 
For the high-temperature expansion also, the presence of the regular term in (j 1 • 2 j) 
makes it difficult to use standard analyzing methods of the series (such as the Pade 
approximation and the inhomogeneous differential approximation) to pick out the 
singular term. However, we believe that the contribution of the singular term to the 
coefficients in the high-temperature series will dominate at sufficiently high orders. 
If this is indeed the case, we could directly compare the high-temperature series of 
the free energy with an expansion series of a form like that of the first term on the 
right-hand side of Eq. (|1-2|) if we calculate a sufficiently long series for the free energy. 
In this paper, we extend the previously derived high-temperature series for the free 
energy—' to a series that is more than two times longer. This makes it possible to 
reveal the singular behavior of Eq. (jl-2p . From an analysis of the obtained series 
we obtain the value of (5 C = 1.1176(33) for the inverse critical temperature and 
b = 1.54(4) for the non-universal constant in Eq. (|l-2p . 

A commonly used method for series expansions is the graphical method. - 1 ' How- 
ever, in this method, one must list all the graphs that contribute to the desired order 
of the series. An alternative, powerful method to generate the expansion series is 
the finite lattice methodP^^'^ It does not require that we list the graphs and 
it reduces the problem to that of calculating the partition functions for the relevant 
finite-size lattices. This is a rather straightforward procedure if we use the transfer 
matrix formulation. In many cases, the finite lattice method generates longer series 
than the graphical methodP^^'^'^'^^® Unfortunately, the original fi- 
nite lattice method can generate a high-temperature series that is at most as long as 
that which can be obtained with the graphical method in the case of the XY model 
in two dimensions. 

Here we apply an improved algorithm of the finite lattice method developed by 
the author and Tabata^S'^ to generate a long series for the the XY model in two 
dimensions. This improved algorithm is powerful in the case of models in which 
the spin variable at each site takes more than two values, including the case that 
it takes an infinite number of values. This algorithm was applied to generate a 
low-temperature series for the absolute value solid-on-solid (ASOS) model and high- 
and low-temperature series for the (/-state Potts model in two dimensions. In both 
cases, it generates much longer series than the original finite lattice method. The 
XY model in two dimensions can be mapped to a kind of solid-on-solid model, and 
the improved algorithm of the finite lattice method enables us to obtain a series 
that is two times longer than the series previously derived for the free energy. We 
confirm from analysis of the obtained long series that the free energy of the two- 
dimensional XY model behaves like Eq. (|1-2|) , with values of the critical temperature 
and the non-universal constant b that are close to the values obtained in studies of 
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the correlation length. 

§2. Algorithm 

We consider the XY model defined on the square lattice. The Hamiltonian of 
this system is 

H = -Y J Wj, (2-1) 

where Sj is a two-dimensional unit vector located at the lattice site i, and the summa- 
tion is taken over all the pairs (i, j) of nearest neighbor sites. The partition function 
at temperature T is 




where 6i is the angle variable of the spin si = (cos0j,sin0j). This model can be 
mapped exactly to a solid-on-solid model whose partition function is 

{h | — oo</ii<+oo} 

where /3 = the quantities are modified Bessel functions, the product is taken 
with respect to all the pairs of neighboring plaquettes, and the variable hi at each 
plaquette i takes an integer value between — oo and +oo. 

The improved algorithm of the finite lattice method to generate the high-temperature 
expansion series for the free energy of this model employed here is the following. We 
first calculate the partition function for each of the l x x l y finite-size rectangular 
lattices with a restricted range of the values of the plaquette variables {h}: 

Z(l x ,ly;h + ,h.)= II^-M^)- ( 2 - 3 ) 

{h\h-<hi<h+}(i,j) 

Here, each plaquette variable hi is restricted as h- < hi < h + (where h- < and 
h+ > 0). We define the size of the finite lattice so that an l x x l y lattice involves 
l x x l y plaquettes, including the bonds and sites on their boundary. For instance, the 
lxl lattice consists of a single plaquette, including 4 bonds and 4 sites. We take 
into account finite-size lattices with l x = and/or l y = 0. An l x x lattice consists 
of l x bonds and l x + ^- sites with no plaquette. The 0x0 lattice consists only of one 
site with no bond or plaquette. We take the boundary condition such that all the 
plaquette variables outside the l x x l y lattice are fixed to zero. 

We then define W of the l x x l y lattice and of the restricted range of the plaquette 
variables recursively as 

W(l x ,l y ;h+,h-) = - log Z(l x ,l y ;h+,h-) 

E V* - l '* + ~ l V + l ) W ^ l V> k '+> ti~)- 

0<l' x <l x , 0<l' y <l y , 0<h' + <h + , h-<h'_<0, 

(2-4) 
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The initial form of W(l x ,l y ; h+,hJ) is as follows: 

Z(0,0;h + ,h^) = 1 and W(0, 0; = 0. 

Z(Z a . l O;/i+,^_) = ^ J 0;0,0)=Jo(^) 1 - for ^ > 1 and 




(2-5) 



(2-6) 



Z(l,l;fc+,/i_) = ^(Pf and 



h=h- 



W(l, 1; 0,0)=- log Z(l, 1; 0, 0) - 2W(1, 0; 0, 0) - 2W(°» 1; 0, 0) - 4W(0, 0; 0, 0) 

= -log/o(/3) 4 + 41og/ (/?) =0. 
W(l, 1; 1,0)=- log Z(l, 1; 1, 0) - W(l, 1; 0, 0) - 2W(1, 0; 1, 0) - 2W(1, 0; 0, 0) 
-2W(0, 1; 1, 0) - 2^(0, 1; 0, 0) - AW (0, 0; 1, 0) - 4W(0, 0; 0, 0) 

= - log{/ (/?) 4 + h (P) 4 } + 4 log J (/3) 



The free energy density / s per site in the infinite-volume limit is given by 



In the standard (graphical) cluster expansion of the free energy for this model, 
a cluster is composed of polymers, and each of the polymers consists of connected 
plaquettes. An example of a polymer can be seen in Fig. 1. A value hi 0) is 
assigned to each site i of the polymer. We can assign to each cluster two numbers, 
Vax (> 0) and h m \ n (< 0), which are the maximum and the minimum, respectively, 
of the quantities hi in all the sites of the polymers of which the cluster consists. 
Then, we can prove^ that the Taylor expansion of W(l x ,l y ; h+,h~) with respect to 




(2-7) 




(2-8) 
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Fig. 1. An example of a polymer in the standard cluster expansion. 



High- Temperature Expansion of the Free Energy in 2D XY Model 



5 



[3 includes the contributions from all the clusters of polymers in the standard cluster 
expansion for which /i max = h + and h m { n = h- and that can be embedded into the 
l x x l y lattice but cannot be embedded into any of its rectangular sub-lattices. The 
series expansion of W(l x , l y ; h+, h~) begins at order f3 n ( lx > l y> h +> h -\ with 

n , h 7 \ f 2{l x + l y ) + 4(h+ + \h-\) - 6 for h+ > and h- < 0, 

n{l x ,l y ,n+,n-)-^ 2 {l x + l y ) + 4(h + + \h-\)-4 for h+ = or h- = 0. 

(2-9) 

An example of a cluster consisting of a single polymer that contributes to the lowest- 
order term of W(l x ,l y ;h+,h-.) is given in Fig. 2. Hence, to obtain the expansion 
series to order (3 , we have only to take into account all combinations of the rectan- 
gular lattice size (l x , l y ) and the range of the plaquette variables (h + , hJ) that satisfy 
the relation n(l x , l y , h + , h-) < N in the summation of Eq. (|2-8|) and to evaluate each 
of the W(l x ,l y ,h+,h-) to order j3 N . 



l_ 


L 


1 




l_ 


L 


h+ 








1 






-1 


1 


1 



Fig. 2. An example of a cluster consisting of a single polymer that contributes to the lowest-order 
term of W(l x , l y ; h+, h-) with l x = l y = 4, h- < and h+ > 0. 

Most of the CPU time required with this algorithm is used for the calculation 
of the relevant partition functions Z(l x ,l y ; h+, hJ) to order (3 N . These partition 
functions are calculated with the transfer matrix method using a procedure in which 
a finite-size lattice is built one plaquette at a timeP^'^ The necessary CPU time 
and memory are proportional to Nl x l y (h + + \h-\ + l) mm fe''a) and N (h + + \h-\ + 
^mm(i x ,i y )^ reS p ec tively. 

If we had applied the original finite lattice method straightforwardly, without 
introducing h + and h- , we would have had to calculate the partition function while 
taking into account all the configurations with each of the plaquette variables ranging 
from —N/4 to N/4 for any necessary size of the lattice. However, because we have 
introduced h + and h-, we only need to take into account the configurations with 
each of the plaquette variable ranging between /i_ and h + , with \h—\ + h + < (N — 
2(l x + l y ))/4. This substantially reduces the number of the configurations that have 
to be taken into account, especially for larger lattices. In fact, if we had used the 
original method, we would have been able to generate a high-temperature series only 
to order /3 24 , instead of f3 . 
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§3. Series 

We have calculated the high-temperature expansion series for the free energy of 
the XY model on the square lattice to order /3 48 . The obtained coefficients for the 
free energy fb per bond are listed in Table I, where the quantities a n are defined 
through the relation 




We have confirmed that each W(l x ,l y ;h + ,h-) in Eq. (|2-8j) starts from the correct 
order in /3, as given by Eq. (|2-9p . The first eleven terms of the series (to order [3 22 ) 
coincide with those obtained by Campostrini et al. P and we have generated thirteen 
new terms. 

To obtain the long numerators and denominators appearing in Table I precisely, 
we have used the following technique. The Taylor series of the Bessel function I n (P) 



Table I. High-temperature series coefficients for the free energy per bond of the XY model on the 
square lattice. 



-1 
-3/4 
-1/9 
31 / 192 
731 / 600 
29239 / 8640 
265427 / 35280 
75180487 / 5160960 
6506950039 / 235146240 
1102473407093 / 26127360000 
6986191770643 / 158070528000 
-1657033646428733 / 49662885888000 
-132067307951625029 / 349709488128000 
-123922162516396049161 / 73815602724864000 
-3988676736069063737749 / 684416646144000000 
-171797935005592667747929271 / 9718351353033523200000 
-4447072728482029389525377383 / 87768860657084006400000 
-1419869374002650746096436008343 / 10036627359845371084800000 
-37739166891866120573923814513909 / 100645068802893860044800000 
-79972307856390845424075199180330349 / 84753742149805355827200000000 
-3529608673175063035175810323068581713 / 1509431550094898846760960000000 
-55921787617500176735150797125026274875543 / 9949597756682679903049482240000000 
-1470445398138160496551515879389714706601279 

/ 119621300301934947016208547840000000 
-7033812372580933320904163849616698204389178239 

/ 287590409630251949139489907015680000000 
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can be written 

^)=£if(§y > (3-2) 

where the coefficients b n j are all integers. The product of two Bessel functions can 
also be written as 

/ ni (/?)/ n2 (/?)=^l(iy , (3.3) 

with 

- i\ 
i'=0 ' V '' 

We note that the combinatorial factor ijx&z%x\ is an integer and that Cj is also an 
integer. Thus, if we treat coefficients like b„ ; i and c, in Eqs. (|3-2p and (|3-3p instead of 
the coefficients of the terms in the Taylor expansion, then we need only treat integer 
values throughout the calculation. For the products and summations of the long 
integers that arise in these calculations, we used the Chinese theorem of modulus.^ 
The calculations were carried out on a workstation at the Information Processing 
Center at OPCT and on an Altix3700 BX2 at YITP of Kyoto University. 



§4. Series analysis 



If the phase transition of the model is of K-T type, the free energy is expected 
to behave as 



/(/3) = A3 2 exp 



2b 



+ B(J3), (4-1) 



(l-/3//3 c )V2_ 

where /3 C is the inverse critical temperature and b (> 0) is the non-universal constant 
that appears in Eqs. f)l • 1 j) and (|l-2p . 

First, it is to be noted that the coefficients of the free energy in Table I exhibit 
repeated changes of sign. This is one indication that the phase transition in this 
model is of the Kosterlitz-Thouless type. This property of the signs of the coefficients 
in the expansion series of the free energy for a model that undergoes a K-T type 
phase transition was first pointed out by Hasenbusch et alP9 The coefficients in the 
Taylor expansion of the first term on the right-hand side of Eq. (|4-ip with respect 
to (3 exhibits repeated changes of sign in general, with the order of the sign change 
depending on the values of c and (3 C , while the coefficients of the terms in the Taylor 
expansion for a quantity with a power law singularity do not exhibit changes of 
sign. In fact, changes in sign of the coefficients is observed in the low-temperature 
series for the surface free energy of the Ising model (i.e., the strong coupling series 
for the string tension of the Z2 lattice gauge theory) in three dimensions^ and 
also in the low-temperature series of the free energy for the ASOS model in two 
dimensions,^ both of which undergo a roughening phase transition of the K-T 
type. 

It is expected that the coefficients of the terms in the series expansion of the 
free energy at sufficiently high orders will be dominated by the contribution from the 
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Taylor expansion of the first singular term in Eq. (14-lh . With this in mind, rewriting 
the high-temperature series as 



N/2 
f = 2 



n=l 




we fit the K-T type function 



f (K ~ T \x) =^xexp 



C + C\X 

fl - xY 



2^ a n 



~(K-T)„TJ 



X 



(4-2) 



(4-3) 



to the high-temperature series under the condition that 



-n 



7, ^2 



E 



N/2 
n=no ' 



(4-4) 



be minimized. In the present case, we have N/2 = 24, and we have carried out the 
fitting for no ranging from 7 to 13. 

We have performed three types of fittings. We first fixed the exponent a to 1/2. 
In this case, the best fitting parameter values are (3 C = 1.1176(33), Co = 2.86(25) and 
c\ = 1.50(18). The uncertainty on each value is due to the range of the values of uq 
adopted. This best fitting value of the inverse critical temperature is consistent with 
the most precise value f3 c = 1.1199(1), obtained from the large-scale Monte Carlo 
simulation,® and the value (3 C = 1.116(4) obtained from the high-temperature series 
for the correlation length and the magnetic susceptibility. 8 '' 7) Although the values 
of Co and c\ have rather large uncertainties, the combination Co + c%, which we find 
to be 4.36(11), has a smaller uncertainty. This is because the value of the quantity 
in Eq. (|4-4p is stable on any line defined by cq + c\ = [constant] near the minimum 
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point. In Fig. 3 we plot a n and a { n ' T) corresponding to the best fitting values of (3 C , 
Co and c\ for no = 13. We find that the agreement between the two series is good not 
only in the range 13 < n < 24 used in the fitting but also in the range 5 < n < 12, 
which was not used in the fitting. This confirms the validity of our assumption that 
the series a n for sufficiently large n is dominated by the Taylor expansion of the first 
term in Eq. (]4- lh . 

Next, we fixed the inverse critical temperature (3 C to 1.1199. In this case, the 
fitting yields a = 0.513(19), and cq + c\ = 4.18(23). This best fitting value of a is 
also consistent with the K-T value of 1/2. 

Finally, we set a to 1/2 and (5 C to 1.1199. In this case, we obtained Co + c\ = 
4.36(11). Here, the behavior of the coefficients a n and a^ < ~ T) with respect to n is 
almost the same as that plotted in Fig. 3. The fitted values of the parameters cq and 
ci imply that the non-universal constant b in Eq. f|l-2j) is given by (c + ci)/2^2 = 
1.54(4). We can see that the value of b obtained from the high-temperature series for 
the free energy is consistent with the value obtained from the Monte Carlo simulation 
of the correlation length to a precision of 15% and with the value obtained from the 
high-temperature series of the correlation length to a precision of 10%. 

§5. Summary 

We have calculated the high-temperature series for the free energy of the XY 
model on the square lattice to order /3 48 using an improved finite-lattice method. 
This method in general enables us to generate a longer expansion series for spin 
systems whose spin variable takes more than two values. The length of the series 
obtained here is two times longer than the series obtained previously with the graph- 
ical method. The results of the analysis for the obtained series give strong support 
to the validity of the prediction given in Eq. ()12j) for the behavior of the free energy 
of the XY model in two dimensions. The value of the critical point derived from 
the series for the free energy is consistent with the value obtained in studies of the 
correlation length employing the numerical simulations and series expansions. The 
value of the non-universal constant obtained from the series for the free energy is 
close to the value obtained from the correlation length. These results reconfirm that 
the phase transition of the model is of the Kosterlitz-Thouless type, which was first 
found in a study of the correlation length. 
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